Resonantly Forced Inhomogeneous Reaction-Diffusion Systems 
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The dynamics of spatiotemporal patterns in oscillatory reaction-difTusion systems subject to pe- 
riodic forcing with a spatially random forcing amplitude field are investigated. Quenched disorder is 
studied using the resonantly forced complex Ginzburg-Landau equation in the 3:1 resonance regime. 
Front roughening and spontaneous nucleation of target patterns are observed and characterized. 
Time dependent spatially varying forcing fields are studied in the 3:1 forced FitzHugh-Nagumo sys- 
tem. The periodic variation of the spatially random forcing amplitude breaks the symmetry among 
the three quasi-homogeneous states of the system, making the three types of fronts separating phases 
inequivalent. The resulting inequality in the front velocities leads to the formation of "compound 
fronts" with velocities lying between those of the individual component fronts, and "pulses" which 
are analogous structures arising from the combination of three fronts. Spiral wave dynamics is 
studied in systems with compound fronts. 



Spatially inhomogeneous reaction-diffusion 
equations may be used to model pattern forma- 
tion and wave propagation in a variety of physical 
systems. We consider situations where the spatial 
inhomogeneity is externally imposed; for exam- 
ple, by inhomogeneous illumination of a light 
sensitive chemical reaction. The focus of the 
investigation is on resonantly forced oscillatory 
systems where the resonant forcing is a applied 
to the system in a spatially inhomogeneous fash- 
ion. Earlier experimental investigations of the 
light sensitive Belousov-Zhabotinsky under spa- 
tially homogeneous resonant forcing by an ex- 
ternal light source revealed phase locked spatial 
patterns. We study wave propagation, pattern 
formation and spiral wave dynamics in oscilla- 
tory reaction-diffusion systems where the applied 
light field has a spatially random intensity pattern 
but varies periodically in time. The phenomena 
we observe include: roughening of fronts sepa- 
rating phase locked domains, nucleation of phase 
locked target patterns and compound fronts with 
distinct properties that give rise to unusual spiral 
wave dynamics. It should be possible to verify the 
phenomena described here by suitably designed 
experiments on light sensitive reacting systems. 



I. INTRODUCTION 

Patterns in reaction-diffusion systems where the kinet- 
ics are spatiotemporally modulated can display a vari- 
ety of pthenomena that are not found in homogeneous 
systemsJd Many reaction-diffusion processes of practical 
interest take place in inhomogeneous media or may be 
coupled to external processes that affect the kinetics in 
a non- uniform manner. A convenient experimental sys- 
tem for studying the effect of spatiotemporal modulations 
on pattern dynamics in spatially distributed systems is 



the ruthenium-catalyzed Belousov-Zhabotinsky (BZ) re- 
actionE This reaction is light-sensitive and thus the ki- 
netics may be modulated by projecting a pattern of illu- 
mination of varying intensity onto the reaction medium. 

Recent studies have made use of the light-sensitive BZ 
system to investigate wavefront propagation in systems 
with spatially disordered excitability. Kadar et al. stud- 
ied stochastic resonance ia a system with a periodically 
regenerated noise pattern.u Sendiha-Nadal et al. studied 
percolation and roughening of wavefronts in sjstems with 
quenched spatially disordered excitability;Q'El excitable 
spiral waves rin the presence time- varying disorder were 
also studied The dynamics of reaction-diffusion waves 
in inhomogeneous excitable media is thought to be rele- 
vant to cardiac fibrillation since electrical waves may bg 
disrupted by irregularities in the heart muscle medium. Q 
Noise is thought to play a role in initiation and propa- 
gation of waves in neural tissue.El Earlier investigations 
into spatially inhomogeneous excitable systems include a 
study of a BZ medium containing catalyst-coated resin, 
beads which served as nucleation sites for wavefronts, u 
numerical simulations of a spatially-distributed network 
of coupled excitable elements which exhibited sponta- 
neous wave initiation—sixuphastic resonance and fragmen- 
tation of wavefronts JIj^EJ and an excitable cellular au- 
tomaton in which the refractory times of the elements 
were assigned randomlyJlj The effect of stochastic spa- 
tial inhomogeneities on other types of reaction-diffusion 
systems has been much less studied. 

Periodically forced reaction-diffusion systems have also 
been investigated. Petrov et al. and Lin et al. subjected 
an oscillatory version of the light-sensitive-BZ-jeaction 
to periodic spatially uniform illumination.E£rE£l As the 
ratio of the forcing frequency to the natural frequency 
neared various resonances patterns were observed. In 
subsequent numerical studies the observed transition be- 
tween labyrinthine and non-labyrinthine two-phased pat- 
terns at the 2:1 resonanfie-iwas reproduced in the periodi- 
cally forced BrusselatorJlxila Belmonte et al. observed a 
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transition from a stable spiral to turbulence in a BZ re- 
action when resonant forcing was applied. O Resonantly 
forced oscillatory systems have been investigated by nu- 
merical sipu-ilation as well as tiieoretically by CcwJlat and 
Emilsson,ll§y CouUet et al.^ Elphick et aZ.,BEl and 
Chate et alx^ These systems exhibit a number of inter- 
esting pattern-forming phenomena. 

Given the wide range of phenomena arising from spa- 
tial disorder in excitable systems and resonant forcing of 
oscillatory systems, one might expect that oscillatory sys- 
tems with spatially inhomogeneous forcing have the po- 
tential to exhibit interesting new features. The research 
presented here explores qualitatively the phenomenology 
of such systems and characterizes some of the phenom- 
ena quantitatively. We restrict the focus of our study to 
systems in two spatial dimensions. 



II. PERIODICALLY FORCED OSCILLATORY 
SYSTEMS 

Consider an externally forced oscillatory reacting sys- 
tem described by the ordinary differential equation, 

^=R(c(t);a,b(0), (1) 

where c{t) is a vector containing the concentrations of 
reagents. The reaction rates are described by the nonlin- 
ear vector function R which depends on a collection of 
parameters a, such as rate constants and constant con- 
centrations of pool chemicals, as well as parameters b(t) 
which comprise the periodic forcing and are of the form 
b(i) = T7g$((jjft) with r7g the constant forcing amplitude 
and $ a 27r-periodic function giving the form of the forc- 
ing. As mentioned above, such forcing may be imple- 
mented by periodic illumination of the system if the re- 
action is light sensitive. 

If h{t) — we suppose the unforced reacting system 
has a stable limit cycle c^{t) with period Tq — I-kIloq. 
In such a system there exists an infinite number of limit 
cycle solutions, c'^ify = Co(t-l- At) which differ from Co(t) 
only by an arbitrary phase shift 2i:lS.t/TQ. Limit cycle at- 
tractors are neutrally stable to phase perturbations cor- 
responding to translations along the orbit. A system fol- 
lowing a limit cycle will, in general, have undergone a 
phase shift when it returns to the limit cycle after expe- 
riencing a small perturbation. 

These characteristics of the unforced oscillaiar may be 
contrasted with those of the forced oscillator .a If LOf/uJo 
is sufficiently close to an irreducible ratio of integers n/m, 
and if the forcing amplitude is sufficiently large, then the 
oscillations may become entrained to the external forc- 
ing and the system possesses n stable limit cycle solu- 
tions of period T = nTf = IwRjiOi ~ mTo which are 
mapped into each other under phase shifts t —^ t + kT/n 
for fc = 0, 1, 2, . . .. A system following one of these limit 
cycles will return to it with no phase shift after a small 



perturbation. This discrete, finite collection of limit cy- 
cles may be contrasted with the infinite and continuous 
collection of limit cycles in the unforced case. The en- 
trained resonantly forced oscillator is a system with n sta- 
ble states, defined by the phase of the oscillations rather 
than by the system's location in phase space. 

The general form of an oscillatory reaction-diffusion 
system with spatially inhomogeneous periodic forcing is 

= R(c(r, t); a, b(r, t)) + BW'c{r, t) , (2) 

where D is a diagonal matrix of diffusion coefficients. 
The parameters responsible for the periodic forcing 
b(r, t) now depend on space as well as time and are of 
the form h{r,t) = rj(r, i)$(wft). The random variable 
ri{r, t) accounts for the fact that the forcing amplitude 
may vary stochastically in space and time. 

In a spatially distributed system with spatially uniform 
forcing, b(r, = h{t) = Tjg$(u;ft), the diffusive couphng 
and the stability of the n limit cycles to phase perturba- 
tions leads to the formation of domains of the different 
phase locked states. At the domain walls separating the 
different phase locked states the phase of the oscillations 
shifts by an amount determined by the character of the 
phase locking. In two dimensions, n-armed spiral waves 
may form, given suitable initial conditions, if the domain 
walls have non-zero velocity. The core of the spiral wave 
is a phase singularity at which all n states meet and 
around which the phase advances by 2m7r. The rotation 
of these phase locked spirals is a result of the propagation 
of the phase dislocations comprising the domain walls; 
thus, they rotate much more slowly than spiral waves in 
the unforced system, where the spiral rotation frequency 
is equal to the frequency of the local oscillations. 

In this article we investigate systems subjected to 
periodic forcing with a stochastic component, either 
quenched disorder, rj(r, t) = rj(r), where the periodically 
applied perturbation has a random distribution in space 
which does not change in the course of the evolution, and 
systems with dynamic disorder, ?7(r, t), where the spatial 
distribution of the perturbation may change with time. 
The reduction of a resonantly forced oscillatory system 
to a corpplex Ginzburg-Landau-iCGL) equation by Gam- 
abaudoE^ and by Elphick et alcB may be simply extended 
to systems with quenched disorder and consequently we 
have employed the CGL equations in our studies of this 
case. For time- varying disorder, one must reconsider the 
derivation of the appropriate CGL equation because of 
the presence of another time scale associated with the 
noise process; hence, we have chosen instead to study 
the FitzHugh-Nagumo system as a typical example of an 
oscillatory reaction-diffusion system. 

In such cases a new length, £s, related to the spa- 
tial correlation range of the noise distribution enters the 
problem. The behavior of the system will depend on the 
magnitude of this length relative to that of other impor- 
tant lengths in the system, such as the diffusion length 
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and the typical width of a domain wall or propagat- 
ing front, Wd- If (-s is sufficiently small compared to both 
of these lengths the system will appear effectively ho- 
mogeneous and behave as if it were subject to periodic 
forcing with amplitude determined by the spatial aver- 
age 77 of r){T). If Is is very large compared to these 
lengths then the system dynamics may be simply rep- 
resented in terms of the dynamics of a collection of uni- 
formly forced patches. The interesting regime is when 
these length scales are comparable and our studies focus 
on these cases. 



For some values of n and 7 Eq. permits multiple i?o 
values, some of which may correspond to unstable fixed 
points. Figure |l| shows a 7 vs. -Ro curve for n = 3; 
the upper branch corresponds to the stable fixed points 
while the lower branch describes the nonzero unstable 
fixed points. The parameter values used to construct 
this figure are jj. — 1, v = \(3\ — 0.6. Note the pres- 
ence of a critical forcing amplitude 7c ~ 0.58 below which 
phase locking does not occur. For n = 3, Eq. (|) gives 
7c = [2((1 + + ^2)) V2 _ + /?^)] 



III. FORCED CGL WITH QUENCHED 
DISORDER 



In the vicinity of the Hopf bifurcation point a periodi- 
cally forced oscillatory reaction-diffusion system may be 
reduced to its normal form,_thfi_forced complex Ginzburg- 
Landau equation (FCGL).OE3 This reduction is valid 
provided the system is near an n : m resonance, where 
n = 1,2,3,4. Such a reduction may also be carried out 
for the case of quenched disorder, ?7(r, t) = T]{r), and the 
FCGL takes the form, 



dA{r,t) 



= (p + iu)A- {l + tl3)\A\^A 



+ 7(r)yl"-^ + (l + ia)VM, (3) 



where 7(r) accounts for the spatial variation of the forc- 
ing amplitude. The complex amplitude A{r, t) describes 
slow modulations in the frequency and amplitude of os- 
cillations of the original system in Poincare planes taken 
at the forcing period; A denotes the complex conjugate 
of A. Equation (|) with 7(r) = 7o, a constant, has been 
used as a model for an oscillatory reaction-diffusion sys- 
tem with spatially uniform resonant forcing.ETtJ 

Consider the normal form of a single resonantly forced 
oscillator obtained by omitting the diffusion terms and 
spatial dependence in 7(r), 



dAit) 
dt 



= {p + iu)A - (1 + i(3)\A\''A + 7A" 



(4) 



For Eq. a critical 7c exists such that for 7 < 7c the 
equation exhibits a stable limit cycle solution, while for 
7 > 7c there are n stable fixed points. These correspond 
to the n stable limit cycle solutions of the original sys- 
tem and can be mapped into each other by phase shifts 
A Ae*^^/". More specifically, the fixed points can 
be found from the stationary solutions of Eq. by ex- 
pressing the complex amplitude in the form A — i?e*'^. 
The modulus, Rq, of the non-zero fixed points of Eq. (Q) 
depends on 7 according to 



+ ml? 



,2(n-2) 




(5) 
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FIG. 1. The modulus, i?o, of the stable (solid line) and 
unstable (dashed line) fixed points of Eq. ^ as a function of 
forcing intensity 7. The parameter values are n = 3, |/9| = 0.6, 
H = 1, V = Q. The three arrows indicate the values 7 = 0.0, 
0.6 and 1.0, which were used in simulations described in Sec- 
The vertical dotted line is located at 7c, the critical 
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tion 

value for phase locking. 

As an example of quenched disorder in resonantly 
forced oscillatory systems the 7(r) fields were taken to 
be dichotomous random variables. Two values for the 
forcing amplitude, 71 and 72, were chosen. The two- 
dimensional system was partitioned into square cells and 
the value of 7(r) in each cell was chosen to be either 71, 
with probability p, or 72 with probability q = 1— p. More 
precisely, if the noise cells have dimension s x s and the 
system's dimensions are W x L — sNw x sNl then 



where 



and 



71 with probability p 

72 with probability q - 



l-p 



(6) 



(7) 
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ey(r) = e.,j{x,y) 

^ e{x- (i - l)s) e{is - x) 
X Oiy - (j - 0ijs 



y), 



(8) 



where 9 is the Heaviside function and (ij) are the discrete 
coordinates of a noise celL For quenched disorder this 
distribution is fixed for all time. The cell size, the values 
of 71 and 72 and the seeding probabilities p and q are 
the relevant parameters to consider. The probabilities p 
and q were typically, but not necessarily, independent of 
position; exceptions will be noted as they occur. This 
7(r) field has the mean value 7 = P7i + 972, and spatial 
autocorrelation 



C(r') = 



(,57(r' + r") Sjjr")) 

1 _ l£l^ (^1 _ if < s and \y'\ < s 

otherwise 



where ^7(r) = 7(r) 
is taken over all r". 



7, r' — {x',y') and the average 



The studies described in this paper investigate pat- 
terns at the 3:1 resonance. For such systems there are 
two interesting cases for resonant forcing with a dichoto- 
mous 7(r) field. In the first case both 71 and 72 lie above 
the phase locking threshold, 7c, so that all regions of 
the medium are entrained to the forcing. In the second 
case only one of 71 or 72 lies above the threshold and 
the medium consists of a mixture of entrained and non- 
entrained regions. In the simulations described below, 




FIG. 2. Interfaces of the 3:1 inhomogeneously forced CGL 
from a single realization in a moving frame, at three well sepa- 
rated times. The phase, cj> of the complex amplitude A — _Re"^ 
is shown, using a gray-scale in which = — tt = +-k is white 
and (p = Ois black. The system size, L x is 800 x 100. The 
noise grain size is s x s = W/25 x W/25. Other parameters 
are given in the text. Boundary conditions are periodic along 
X and no-flux along y. 



Front roughening is also observed if 71 < 7c < 72i (i-e. 
(9)when the medium consists of a mixture of entrained and 
non-entrained regions) but is difficult to study because 
spontaneously nucleated patterns interfere with propa- 
gating fronts. This case will be examined below. 

The propagating fronts in this system experience lo- 
cal velocity fluctuations arising from spatial variations 
in the 7(r) field. Diffusion will tend to eliminate front 
roughness generated by random fluctuations in the 7(r) 
field; consequently, the front dynamics, should obey the 
Kardar-Parisi-Zhang (KPZ) equation,cZl 



dh{x, t) 
dt 



D 



d^h \ fdh 



dx^ 2 \dx 



+ C(a:,t) 



(10) 



including those in Sec. IV, numerical integration was per- 
formed using explicit forward differencing and a second- 
order discrete Laplacian. 



A. All sites phase locked; front roughening 



When both 71 and 72 lie above the phase locking 
threshold 7c then all regions of the medium are tristable. 
The system possesses three quasi-homogeneous station- 
ary states in which the complex amplitude fluctuates 
about an average value. Domain walls separating these 
phase locked states are in general non-stationary in the 
3 : 1 resonant regime since all three phases are inequiv- 
alent, although their velocity may pass through zero as 
parameters are tuned. Initially planar fronts in these 
inhomogeneously forced systems roughen as they prop- 
agate. Figure || shows an example of a rough interface 
separating two of the three phases. 



where h{x,t) is the front position, v is the average front 
velocity, D and A are phenomenological coefficients and 
^{x, t) is Gaussian white noise with zero mean and cor- 
relation {C{x',t')(:{x",t")) = 2r S{x' ~ x")S{t' - t"). 
In such a circumstance the interface width w(t) — 

{L^^ J2ii^i^iT^) ~ ^(0)^ increases with time as 

w{t) ~ for short times; the average width of a satu- 
rated front, Ws, scales with system size as Ws L" , where 
a = 1/2 and (3 = 1/3. 

We have verified these KPZ scaling properties for a 
FCGL system with parameter values {a = 1, (3 = 0.6, 
fj, + iv = 1) and the forcing field parameters (71 — 0.60, 
72 = 1, p = 0.50) for which the critical forcing ampli- 
tude is 7c « 0.58 (cf. Fig. |l|). The front properties were 
measured in a frame moving with the front. The system 
dimensions were W = 100 and L = 100, 200, 400 and 
800. The noise grain size was sx s = W/25 x W/25. The 
boundary conditions were periodic along the boundaries 
perpendicular to the front, and no- flux along the bound- 
aries perpendicular to the direction of front motion. 

The saturated front width Ws was found to scale with 
system size as Wg ~ L^/^. Plotting {w{t))/L" against 
t/L"^^ collapses the {w{t)) versus t data for four differ- 
ent L onto a single curve when the KPZ values a — 1/2, 
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/3 = 1/3 are used, as shown in Fig. 
an average over realizations. 



Here (•) denotes 



The width temporal autocorrelation function 



(Sw{t + t') Sw{t')) 
{dw{t') 6w{t')) 



(11) 



was calculated for fronts in the saturated regime at four 
system sizes L = 100, 200, 400 and 800. Here (•) repre- 
sents a time average within the saturated regime. Tem- 
poral correlations were found to decay to zero, with the 
rate of decay decreasing as system size increases (Fig. ^). 
This is expected, since the larger system size allows larger 
fluctuations in the front profile which require longer times 
to form and decay. 
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FIG. 4. The width temporal autocorrelation function for 
saturated interfaces in the 3:1 inhomogenously forced CGL 
for system sizes L = 100 (dotted line), 200 (short dashes), 
400 (long dashes) and 800 (solid line). 
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To investigate the presence of spatial correlations in 
saturated front profiles, the height spatial autocorrela- 
tion function 



was calculated for the various system sizes. In Eq. ( [T^ ) 
the inner (•) refer to averaging over x' in a single front 
profile and the outer (•) refers to averaging over different 
saturated front profiles. Results for L = 100 and L ^ 800 
are shown (Fig. ^). 



FIG. 3. Early time {w{t)) vs. t curves for initially planar 
fronts in the 3 : 1 inhomogenously forced CGL rescaled ac- 
cording to the KPZ scaling exponents. Curves are shown for 
system sizes L = 100 (dotted line), 200 (short dashes), 400 
(long dashes) and 800 (solid line). Each curve results from 
the average over 80 realizations of the stochastic dynamics. 



In the case of a front free of spatial structure, i.e. a 
periodic random walk that returns to its initial position 
in N steps, the height spatial autocorrelation function is 
of the form Cf (x) = 1 - (2D/L)x(L -^), where D is a 
phenomenological diffusion coefficient. E3 The deviations 
of the fit of the numerically determined Ch{x) from a 
quadratic function indicate the existence of spatial cor- 
relations. 
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FIG. 5. The height spatial autocorrelation functions for 
saturated fronts in systems of size L — 100 (solid line) and 
L = 800 (long dashes). The best fit quadratic functions are 
shown for L = 100 (dotted line) and L — 800 (short dashes). 



B. A medium with phase locked and oscillatory 
sites; spontaneous nucleation of target patterns 



When 7i < 7c and 72 > 7c, the system consists of a 
mixture of tristable regions and osciUatory regions. If the 
density of oscillatory sites, p, is low the diffusive coupling 
maintains the value of A within the oscillatory regions 
near that in the adjacent tristable regions. The medium 
behaves essentially like a tristable medium and supports 
three quasi-homogeneous stable stationary states, travel- 
ling kink-like domain walls and three-armed spiral waves. 
The oscillatory sites provide a spatial inhomogeneity that 
leads to roughening of domain walls and to fluctuations 
of the concentrations within domains. 

With increasing p, target patterns are observed 
(Fig. They consist of concentric, approximately cir- 
cular domain walls moving outwards from a central re- 
gion (a "pacemaker") where they are initiated periodi- 
cally. Within each ring of the target the concentration 
is quasi-homogeneous. Note that all images in Fig. || are 
from realizations with identical parameters, however a 
range of wavelengths is observed. 




FIG. 6. Target patterns generated in the 3:1 forced CGL 
from spatially homogeneous initial conditions. The phase, 
(j) of the complex amplitude A = Re"^ is shown, using a 
gray-scale in which <j!> = — tt = +n is white and <^ = is 
black. Forcing field parameters are (71 = 0, 72 = 1, p = 0.30). 
Other parameters are {fi + iv = 1, a — 1, P — 0.6), for which 
7c ~ 0.58. The system size is L x L = 200 x 200, the noise 
grain size is s x s = L/200 x L/200. Boundary conditions are 
periodic. 

The probability of a realization possessing a target pat- 
tern was measured as a function of p and system size 
(Fig. 1^). It was found to approach zero for low p, one 
for high p, and to increase rapidly around some criti- 
cal pc- Figure |^ shows results for the parameter values 
{fi -\- iv — 1, a ~ 1, P = —0.6) and forcing field param- 
eters (71 = 0, 72 — 1); qualitatively similar results were 
also obtained for {fi + iv = 1, a = 1, P = 0.6, 71 = 0, 
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72 = 1) and {fj, + iiy = Q, a = 2, (3 = -0.82, 71 = 0, 
72 = 1). For larger system sizes, the probability of oc- 
currence of a target pattern is higher, consistent with a 
fixed probability per unit area of the medium nucleating 
a target pattern. 




0.1 0.15 0.2 0.25 0.3 0.35 

P 

FIG. 7. The fraction of realizations in which one or more 
pacemakers occurred as a function of p, the density of oscil- 
latory (7 = 71 = 0) sites. Results from numerical simulations 
with noise grain size s x s — 0.5 x 0.5 are shown for system 
sizes 50 X 50 (a) and 100 x 100 (▼). The curves predicted 
by the model described in the text (Eqs. (|l|)-(|2o|)) with the 
parameter values 7* — 0.65 and 3sfp = 180 are shown for these 
same system sizes as solid and dashed lines, respectively. 

The occurrence of target patterns may be explained 
by supposing that diffusion provides a mechanism for av- 
eraging 7(r) over some length scale, and that regions of 
the medium behave qualitatively like a uniform medium 
with the same local average 7(r). Thus the medium is 
locally tristable for a high 7(r), but if 7(r) is sufficiently 
low it will be oscillatory. As the density of oscillatory 
sites increases, it becomes increasingly likely that there 
will exist regions with a local 7 sufficiently low to exhibit 
oscillatory dynamics. 

We can determine a profile of the average 7 field around 
the average pacemaker. We define 7p as the average value 
of 7 a distance R from the center of a pacemaker. This 
is given by 



d_ 

2TrRdR 



where 



ttR^ 



7(r) dr 



(13) 



(14) 



|r-ro|<-R 



parameters. In practice, 7p is calculated as the discrete 
derivative 

^P(^)^ n{R + ARr-nR^ • 

Figure || shows % and 7d for p = 0.30, system size 
L X L = 100 X 100, noise grain size s x s = L/200 x L/200, 
and other parameter values equal to those used to mea- 
sure the nucleation probability curves shown in Fig. |^. 
Qualitatively similar plots, in which 7p and 7d increase 
from a low value at i? = to the mean-field value, 7, at 
high R were observed for other values of p. Also shown in 
Fig. H are 7p calculated from the neighborhoods of points 
chosen randomly from the realizations rather than from 
points centered on pacemakers, and the mean-field value, 
7- 




is the average value of 7 over a disk of radius R, averaged 
over all pacemakers in all realizations with a given set of 



FIG. 8. The solid curve shows 7p(i?); the long-dashed 
curve is 7d. Averages are taken over 20 realizations. The 
short-dashed curve shows 7p calculated around points chosen 
randomly from the realizations. A horizontal line indicates 
7 — 0.70. The values 7* — 0.65, 3? = 3.8, calculated using the 
pacemaker model, are also shown. 

Consider the following simple model of pacemaker for- 
mation. We divide the system into 1^ "sites" of size 
Ax X Ax such that the value of 7 is constant over a 
site; i.e., the sites may be noise domains or subdivisions 
of noise domains. We assume that, on average, pacemak- 
ers have a radius of Jl and contain JvTp — TT{Ji/Ax)'^ sites. 
Thus we divide the system into J^/'Np independent do- 
mains which are potential pacemakers. We assume that 
such a domain is a pacemaker if the average value of 7 
within it is less than some threshold for oscillatory be- 
havior, 7*. If Ni is the number of sites on which 7 = 71 
and N2 = — Ni is the number of 72 sites, then the 
potential pacemaker is a pacemaker when 



7 



7iA^i + I2N2 _ liN^ + 72(?^p - N2) 



that is, when 
Defining 



72 - 71 



< 



7 , 



N* = 



72 ~ 7 



72 - 71 

the probability that a domain is a pacemaker is 



(16) 



(17) 



(18) 



(19) 



and it is not a pacemaker with probabihty Q = 1 — J". 
The probability that the entire system possesses no pace- 
makers is Q-'^Z-'^p and it possesses one or more pacemakers 
with probability 



P(3^,p) = 1 - Q^/^" 



(20) 



The parameters 3\fp and 7* uniquely determine the 
P(!N', p) surface. The values giving the best fit to the 
experimental data were 7* = 0.65 and 3\fp = 180 which 
implies 3? 3.8. These values are compared with the nu- 
merically determined jd{R) curve in Fig. ^ The simple 
model predicts that the average pacemaker should have 
7d(3?) = 7*; the point (3?, 7*) lies close to, but not on, 
the 7d curve. 

In order to provide further insight into the criteria nec- 
essary for a region to act as a pacemaker, a series of 
studies was carried out in which the 7 field consisted of 
a disk of radius R sites with a density pin of oscillatory 
sites. This disk was embedded in a field with density of 
oscillatory sites Pout- In all cases pi^ was greater than 
Pout, so that the disk could act as a pacemaker. Multiple 
realizations of the evolution were simulated at various 
values of R, pin and Pout, and the fraction of realizations 
Pr(br) in which the central disk emitted target waves 
into the surrounding medium (i.e. in which "breakout" 
occurred) was measured. 

Figure ^ shows Pr(br) as a function of R and pin for 
Pout = 0.10. As one would expect, Pr(br) increases as R 
increases and decreases as pin decreases, with the excep- 
tion that at Pin = 1 the probability of breakout increases 
with R until i? = 3 and then decreases for 3 < i? < 5, 
after which it increases, (cf. Figs. I and 0). For 
Pin = 0.95, when there is a small fraction of tristable 
sites in the inner disk, the magnitude of the decrease 
is substantially reduced (cf. Fig. [lO| ). Apart from the 
anomalous region at low pin the trends that Pr(br) in- 
creases with increasing R and decreases with decreasing 
Pin are consistent with the notion that pacemakers form 
when the local density of oscillatory sites is high. The 
behavior for pout — 0.15 was qualitatively similar. 




FIG. 9. The fraction of realizations in which target waves 
were emitted from the central disk region. Parameter values 
apart from p are as in Figs. |^ and ^ the noise grain size is 
0.5 X 0.5. Here pout = 0.10. 



Pr(br) 




FIG. 10. Fraction of realizations in which target waves were 
emitted from the central disk region as a function of R, for 
Pout = 0.1 and Pin — 1.0 (solid line) and pin = 0.95 (dashed 
line) . 

The anomalous decrease between i? = 3 and i? = 5 is 
possibly related to the fact that the anomalous behavior 
begins when R^i £u, where £d is the diffusion length of 
the unforced system (7(r) = 71 =0). We find the dif- 
fusion length io — V Dt « 3.24 by taking the diffusion 
coefficient D to be unity, and the characteristic time r 
to be equal to 27r//3/x, the period of homogeneous oscilla- 
tions in the unforced system for the parameters /3 = —0.6, 
^ = 1 used in these studies. Parts of the system sepa- 
rated by a length greater than £]j evolve independently 
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over time intervals less than r. For large R one observes 
that pacemaker nucleation occurs locally on the bound- 
ary of the disk. As the disk perimeter increases with R 
so does the probability of forming a local pacemaker on 
the disk boundary. 

In addition to the target patterns discussed above, one 
may also observe spiral waves if the initial conditions 
contain a phase defect. An example of such a spiral is 
shown in Fig. |ll]. It was formed in a realization with 
7 = 0.50 < 7c, hence the medium may be thought of as 
oscillatory, exhibiting three-fold symmetric relaxational 
local dynamics, rather than as tristable. 




is updated on a time scale that is some multiple of the 
forcing period. 

Before considering Eq. ( pT] ) we discuss the associated 
system of ordinary differential equations 

—u(t) — u — u ^ V 
dt ^ ' 

^u(t) = e(M - a-y + 6) , (22) 

in which 6 is a constant. If < a < 1 the system pos- 
sesses a single fixed point. If, in addition, 6 = 0, then 
the fixed point is unstable and the system exhibits a 
stable limit cycle whenever ae < 1. In this article we 
consider only systems with < a < 1 and ae < 1. Fig- 
ure ^ shows the limit cycle for a system with a = 0.3, 
e = 0.1 and & = 0. Also shown are the cubic u-nuUcline 
and the linear ti-nuUcline for different h values. The ef- 
fect of varying h is to shift the i;-nullcline relative to the 
M-nuUcline. As \h\ increases from zero the limit cycle 
contracts in the phase plane, eventually collapsing to a 
stable fixed point at a Hopf bifurcation that occurs when 

|6| = hH^ (l-(2a + a2e)/3)((l-ae)/3)^''^ For I6| > hu 
the system exhibits excitable kinetics. Figure |l^ shows 
w-nuUclines corresponding Xo h — ±0.46, which lies just 
inside the excitable regime. 



FIG. 11. Spiral wave in the 3:1 inhomogeneously forced 
CGL. The phase, of the complex amplitude A — Re^'^ is 1.5 
shown, using a gray-scale in which — — tt = +Tr is white 
and <j) — is black. Parameters: a = 1, f3 — —0.6, /i + iu — 1, 
71 = 0, 72 = 1, p = 0.50, L = ICQ, s = L/200, no-flux 
boundary conditions. 



V - 

IV. SYSTEMS WITH TIME- VARYING SPATIAL 
DISORDER 



We now consider situations where the spatial distribu- 
tion of forcing amplitudes varies in time. For this pur- 
pose we examine the behavior of the spatially distributed 
FitzHugh-Nagumo (FHN) system 



du{r, t) 

dt 
dv{r, t) 

dt 



e{u - av + 6(r, <)) + DyV'^v 



(21) 



subject to such time- varying noise distributions. Here 
M(r,t) and w(r, i) are the "concentrations", a and e are 
constant parameters, Du and Dy are the diffusion co- 
efficients and &(r, t) is a forcing function of the form 
?7(r, t) cos LOit with the forcing amplitude ?7(r, t) a random 
variable. We have chosen this form to make explicit the 
periodic component of the external forcing whose ampli- 
tude is a periodically or randomly updated spatial ran- 
dom variable. In the applications described below rj{T, t) 



■1.5 



■1.5 





u 



1.5 



FIG. 12. The nuUclines and limit cycle of the FHN system 
(Eq. |2^ ), with a = 0.3. The cubic u nuUcline and linear v 
nullcline corresponding to 6 = are shown as solid lines. The 
dashed lines are the v nuUclines for b = ±0.46, for which the 
system is excitable. The limit cycle for 6 = 0, e = 0.1 is also 
shown. 

Given this background, we consider the dynamics of 
Eq. ( |2l] ) as an example of an oscillatory reaction-diffusion 
system with periodic forcing. The forcing amplitude field 
?7(r, t) was similar to the 7(r) fields used in the quenched 
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disorder studies described earlier in that the system was 
divided into squares which were randomly assigned one 
of the two forcing intensities rji and 772 . This disordered 
forcing amplitude field was periodically updated, with 
the new values of the amplitude in the spatial distri- 
bution drawn from the same dichotomous distribution 
of amplitudes. The updating period was taken to be 
nTf time units, which is the period of the corresponding 
n : m entrained ordinary differential equation. (The in- 
vestigations described here consider only the case where 
the updating is on resonance, i.e. where the interval be- 
tween updates is nTf. We shall not consider the case of 
periodic updating where the updating is off-resonance.) 

The nTf-periodic updating of the forcing amplitude 
field may have a phase offset relative to the Tf-periodic 
forcing coswft. We describe this offset with the param- 
eter a which ranges from to 1 and specifies the phase 
offset in units of nTf, i.e. as a fraction of a period. Thus, 
updates occur at times 

Tk ^ {ik-l) + a)nTf for fc = 1, 2, 3 . . . , (23) 



in addition to the initial specification of the random forc- 
ing field at To = 0. To formalize the foregoing: the studies 
were carried out with 6(r, t) = 7/(r, t) cos ujft, where 

00 Nnr Nl 
fc=0 1=1 j=l 

(24) 

where 9 is the Heaviside function, Oij(r) is the charac- 
teristic function selecting the square with discrete coor- 
dinates (ij) as in Eq. (||), and 

^ r ?7i with probability p 

1 ri2 with probability q — 1 ~ p . ^ ' 

The spatial average and time average are equal and given 
by ^{i) = PVi + = ('7(r)). The space-time autocorre- 
lation function is 



^^^^^^^{Svi'r' + r,t' + t) Sv{r',t')) 



(26) 



(1 -¥)(!- ¥)(! - A) if 1^1 < ^ \y\ < ^ and \t\ < uT, 



otherwise 



The investigations described in this section considered 
systems at the n : m = 3 : 1 resonance, using the param- 
eters a = 0.3, e = 0.1, Wf/wo = 3.05, = = 0.25 
with the forcing field parameters 771 = 0, 772 = 0.92, 
p — q = 0.50 and noise grain size s x s = 4x4. 
The corresponding mean field system, Eq. ( |2l| ) with 
6(r,t) = ryocoswfi, rjo = fj = 0.46, lies in the entrained 
regime and admits three-armed phase locked spiral waves 
(Fig. |ll). 

Figure ^ shows an example of a spiral wave in the 
FHN system with quen ched disorder, analogous to that 
considered in Sec. IIIB for the CGL equation. For this 



case the ry(r, t) field may be described within Eqs. (|24|)- 
( p5| ) if we take ri = -l-oo. Substantial front roughening is 
apparent. This system exhibits a phenomenon not seen in 
the studies in the inhomogenously forced CGL described 
in Sec. III. In the uniformly forced FHN, with 77(r, t) = tjq 
and other parameters the same, the front velocity passes 
through zero as 770 is varied. Variations in the effective 
local rj values, combined with front curvature effects, re- 
sult in frequent local pinning of the fronts. The fronts 
may be depinned through coupling to mobile portions of 
the front, or by the perturbation provided by a following 
front as it approaches near the pinned front. Thus, the 
fronts move with an irregular stop-start motion that is 
controlled by the pinning and depinning events. It seems 



unlikely that the resulting fronts could obey KPZ scaling; 
indeed, inspection of Fig. |l^ suggests the front profile is 
unlikely to even remain a single-valued function of po- 
sition. Realizations of spirals in a system with smaller 
size eventually reached a stationary configuration where 
fronts were pinned everywhere along their lengths. 




FIG. 13. A three-armed spiral wave in the forced FHN 
reaction-diffusion system (Eq. (|2^)) with spatially uni- 
form {ri{r) = rjo — 0.46) forcing near the 3:1 resonance 
{uf/cuQ — 3.05). The gray-scale indicates the value of 
tan^^(ii/M). The system size is 512 x 512; boundary con- 
ditions are no-flux. 



10 



For the FHN system with time-varying disor- 
der and the aforementioned parameters, three quasi- 
homogeneous states were observed, similar to the behav- 
ior seen in the CGL with quenched disorder. The ex- 
istence of noise-update events breaks the symmetry be- 
tween the different entrained states of the system. Simi- 
larly, domain walls are now no longer equivalent and may 
travel at different velocities. Consider the 3:1 forced sys- 
tem and arbitrarily label the phases 1, 2, and 3. In the 
following discussion, a [31] front means a domain wall 
between phases 3 and 1, with phase 3 on the left; its 
opposite front is [13]. There are three front types: [31], 
[12], [23] (and their opposites). The velocities of these 
fronts were measured as function of a (Fig. It suf- 
fices to measure the velocities for 5/6 < cr < 1, the val- 
ues for other a follow from the system's symmetry under 
t ^ i + Tf and (u, v, t) {-u, -v, t + Tf/2) . AU three 
fronts move to the left (positive velocity). 




FIG. 14. A three-armed spiral formed in the inhomoge- 
neously forced FitzHugh-Nagumo system with quenched dis- 
order. The parameters are identical to those in Fig. ^ ex- 
cept for the forcing field parameters which are ?7i = and 
772 = 2770 = 0.92 and p = q — 0.50. For this forcing 
field f] — rjo = 0.46. The gray-scale indicates the value of 
taii^^{v/u). The system size is 512 x 512; boundary condi- 
tions are no-flux. 



Depending on a, their velocities rank as V12 > V23 > 
V31 or as V23 > V12 > W31. We note that for all cr, 
1^12 > 1^31- Thus, if a system has initial conditions consist- 
ing of two plane fronts, [31 ... 12] (where . . . represents a 
region of phase 1), the [12] front will move faster than the 
[31] front and the distance between the two will decrease. 
Eventually, the [12] front will closely approach the [31] 
front and a new stable propagating front consisting of a 
thin layer of phase 1 connecting phases 3 and 2 will re- 
sult. We term such a front a compound front and denote 
it [312]. Similarly, for a values where V23 > V12 > W31, 
the compound front [123] exists and can be obtained from 
the starting configuration [12 . . . 23]. 




0.028 



0.024 



0.020 



0.016 



FIG. 15. Front velocity versus a in the forced FHN with 
periodically updated spatial disorder. Lines of the thinnest 
width indicate simple fronts: [23] (long dashes), [31] (short 
dashes), [12] (solid line); medium width lines indicate com- 
pound fronts: [123] (solid line, exists for 0.99... < cr < 1), 
[312] (dashes); thick lines indicate pulses: [3123] (dashes, ex- 
ists for 0.985... < cr < 1), [2312] (dotted line, exists for 
5/6 < a < 0.985...). In addition, for 5/6 < cr < 0.848..., 
where W31 > V23 there should exist a [231] compound front; 
this front has not been characterized. Velocities were mea- 
sured in a moving frame in a 200 x 200 system; the average 
front position was measured stroboscopically with period nT{ 
and linear regression was performed to find the slope. 

As one might expect, compound fronts cannot be made 
from a slower moving front following a faster moving 
front. For example, [231] is not stable; it splits into 
[23 . . . 31] because V23 > W31. 

The velocities of these compound fronts were measured 
as a funtion of a and, within our numerical accuracy, were 
found to lie between the velocities of the two simple fronts 
from which they were derived (i.e. W12 > W312 > W31 and 
V23 > V123 > V12 ) (Fig. |l^). Depending on a either 
V312 > V23 or V23 > V312. For V312 > V23 one expects 
the travelling pulse solution [2312] to be stable, and this 
is indeed the case, while for V23 > W312 the pulse [2312] 
is unstable (it splits into [23] and [312]) while the pulse 
[3123] is stable. The pulse velocities were measured, they 
are essentially the same as the velocity of the faster mov- 
ing component of the pulse. 

The existence of stable pulse solutions joing two do- 
mains of the same phase raises the question of whether a 
"one-armed" spiral whose arms consist of the pulse can 
exist. Figure 16 shows a stable spiral for cr — 11/12, a 
regime where the velocity ordering is f 12 > ^312 > V23 > 
V31. It is not a "one-armed" spiral with a [2312] pulse 
front but could be viewed as a two-armed spiral with 
arms consisting of phases 3 and 2, with fronts of type 
[23] and [312]. Since the [312] front velocity is greater 
than that of the [23] front, one expects that as the waves 
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travel outward phase 3 will shrink and phase 2 will grow, 
and far from the core the waves will become a train of 
[2312] pulses. 




FIG. 16. A spiral wave in the forced FHN with 
time-varying spatial disorder with a — 11/12. The gray-scale 
indicates the value of tan~^(w/u). The phases are 1 (light 
grey ), 2 (dark grey), and 3 (medium grey). The system size 
is 1024 X 1024. The noise grain size is s x s = L/256 x L/256. 
Boundary conditions are no-flux. 

Figure ^ shows a spiral for a = 1. In this regime the 
velocity ordering is V2z > V123 > W3123 > V12 > "312 > 
?;3i. The stable pulse is [3123] rather than [2312]. Far 
from the core we expect the waves to become a train of 
[3123] pulses and in Fig. ^ this can indeed be seen to 
happen. 




The noise grain size is s x s = L/256 x L/256. Boundary con- 
ditions are no-flux. 

The motion of the spiral core was recorded for a realiza- 
tion of the dynamics with a = 11/12. The trajectory of 
the core, rdt) = {xc{t),yc{t)), describes a "noisy flower 
pattern" (Fig. |l^) . Both the periodic looping motion and 
distortions of the simple flower pattern due to the noise 



are evident. A plot of (|rc|^) vs. t shows periodic behav- 
ior with period '--^ 17 000, which is also the mean period 
of rotation of the spiral. In the mean field system with 
i]{r,t) = fj the core is stationary. It is also stationary 
for the uniformly forced systems with r](r, t) = rji and 
T](r,t) = 772, the two extreme values of 77 in the dichoto- 
mous noise process. Consequently, the core motion is a 
result of the time-varying spatial disorder of the forcing 
amplitude field. 




FIG. 18. Space-time plot of spiral core position ver- 
sus time for a realization of the forced FHN dynamics with 
time-varying spatial disorder in a 512 x 512 system. The up- 
dating parameter is cr = 11/12. 

We have also investigated time- varying noise where up- 
dates occurred at Poisson-distributed intervals instead of 
periodically. The Poisson distribution used was Pr{t < 
Ate < t + dt) = (l/t)e-*/* where t = {Ate) = nTf. 
With this choice of t the mean time between updates is 
the same as in the on-resonance periodic updating case 
discussed previously and corresponds to one period of the 
entrained system. Thus, if Ate are chosen from this dis- 
tribution then updates of the forcing field occur at times 

k 

e=i 

in addition to the initial specification of the forcing field 
at tq. With these definitions of Tk in place of Eq. (|2^), 
7j{r,t) is as given in Eqs. (p^-(p5|). 

We expect that in this system the three phases will be 
equivalent on average on time scales longer than the av- 
erage interval between r](T) field updates. The observed 
spiral shown in Fig. |l^ confirms this equivalence. The 
three arms seen in the figure are approximately equiva- 
lent and, when the animation of the dynamics is viewed, 
this equivalence is preserved in time. 



12 




FIG. 19. A spiral wave in the forced FHN system with 
time-varying spatial disorder in which the interval between 
updates is chosen from a Poisson distribution. The gray-scale 
indicates the value of ta.n~^{v/u). The system size is 
1024 X 1024. The noise grain size is s x s = L/256 x L/256. 
Boundary conditions are no-flux. 



V. DISCUSSION 

We have explored the phenomenology of resonantly 
forced oscillatory reaction-diffusion systems subject to 
both quenched and time-varying disorder in the forc- 
ing amplitude field. Noteworthy phenomena found when 
there is quenched disorder are front roughening and spon- 
taneous nucleation of target patterns. Spontaneous nu- 
cleation of target patterns arises because diffusion effec- 
tively causes averaging of 7(r) ocally over some length 
scale; hence, the medium locally behaves like a uniform 
system with the same 7. Alternatively but equivalently, 
we may describe the dynamics as arising from compe- 
tition between two regimes selected by the dichotomous 
forcing values which lie on either side of a bifurcation 
point. In some spatial regions one regime dominates the 
dynamics while in other regions the second regime dom- 
inates. 

For time-varying disorder, the symmetry of the sys- 
tem's three states is broken, and thus the velocities of 
different domain walls differ. As a result, travelling front 
structures other than simple kink-like fronts may exist. 
These are the compound fronts and pulses. The veloci- 
ties of the domain walls (lepc;nd on the phase of the forc- 
ing field updating, therefore the updating parameter a 
selects between regimes in which different sets of com- 
pound fronts and pulses exist. This asymmetry among 
states also leads to spiral waves with inequivalent arms. 

Meandering core motion with a noisy flower-like trajec- 
tory is seen with time- varying disorder. This core motion 
arises purely from inhomogeneities in //(r, t) since none 
of the uniformly forced systems with r]{r,t) = j7,r?i,r?2 
exhibit core meandering. 

Some of the effects described herein may arise from the 
spatial inhomogeneity of the forcing, and not necessarily 



from its disordered nature. One may also consider regu- 
lar patterns of inhomogeneous forcing and investigations 
of this type are in progress. Many of the effects found in 
this study, such as target pattern nucleation, KPZ front 
roughening, and front dynamics for time- varying disorder 
arise from the stochastic nature of the forcing field and 
would not be found in a system with a regular pattern of 
forcing. 

There are opportunities for further exploration of spa- 
tially disordered resonantly forced systems; for example, 
the effects of noise on various phenomena or bifurcations 
known in spatially uniform resonantly forced systems 
have not been studied. Additionally, one may investigate 
oscillatory systems possessing richer limit cycle dynamics 
than the relatively uncomplicated examples considered in 
this paper. 

Systems of the types described here could be readily re- 
alized in an experimental setup. Standard methodology 
for investigation of spatial disorder in the light-sensitive 
excitable BZ reaction involves use of a computer- 
controlled video projector to project a precisely con- 
trollable spatiotemporal pattern of illumination inten- 
sity onto the reaction medium. Thus, the only neces- 
sary changes are the use of the light-sensitive BZ in the 
oscillatory regime and reprogramming of the projector 
to provide a periodic illumination signal incorporating 
appropriate stochastic spatiotemporal modulation of the 
light intensity. 
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